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§ ■ ABSTRACT 
(N 

We examine the growth of eccentricities of a population of particles with initially 
' circular orbits around a central massive body. Successive encounters between pairs of 

OO ' particles increase the eccentricities in the disk on average. As long as the epicyclic 

motions of the particles are small compared to the shearing motion between Keplerian 
orbits, there is no preferred scale for the eccentricities. The simplification due to this 
self-similarity allows us to find an analytic form for the distribution function; full nu- 
\ merical integrations of a disk with 200 planetesimals verify our analytical self-similar 

distribution. The shape of this non-equilibrium profile is identical to the equilibrium 



00 

' profile of a shear-dominated population whose mutual excitations are balanced by dy- 

namical friction or Epstein gas drag. 
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Q^' Subject headings: planets and satellites: formation — solar system: formation 

o 

S : 1. INTRODUCTION 



Modern computational power allows the simultaneous integration of the orbits of increasingly 
numerous particles. Much of the planet formation process, however, involves particle numbers that 
exceed the limits of computational efficiency. This limitation is often circumvented with a statistical 
approach. By monitoring the gravitational interactions of the particles in a time-averaged sense, 
various properties of the particle population can be calculated without a full A^-body simulation. 

Collins & Sari (2006), hereafter Paper I, motivate a Boltzmann equation to describe the evolu- 
tion of the eccentricity distribution of an ensemble of particles in which the relative motion between 
any two interacting particles is dominated by the shearing motion of close circular orbits. Such a 
regime of orbital eccentricities is called shear-dominated. The solution of their equation provides a 
simple analytic expression for the equilibrium eccentricity distribution that results when dynamical 
friction can balance the mutual interactions of the particles; the analytic expression matches results 
from numerical simulations remarkably well. 



In this letter we derive analytically the non-equilibrium distribution function of interacting 
shear-dominated particles in the absence of dynamical friction. §2 reviews the construction of the 
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Boltzmann equation. In §3, we show that the distribution function behaves self-similarly, and the 
shape of the non-equihbrium distribution function is identical to the equilibrium distribution of 
Paper I. §4 corroborates this result with numerical simulations. Conclusions follow in §5. 

2. THE TIME-DEPENDENT BOLTZMANN EQUATION 

We consider a disk of particles on initially circular orbits around a massive central body. We 
write their surface mass density a and the mass of a single body m. The number density, a/m 
is sufficiently low that three-body encounters are very rare, therefore the orbital evolution of each 
body is well described as a sequence of pair-wise encounters. 

The change in eccentricity due to one such encounter can be calculated analytically. For 
completeness, we summarize the derivation presented in Paper I. Let one particle, with a semi- 
major axis a, encounter another with semi-major axis a + h. In the limit of 6 <C a, the relative 
orbital frequency between the pair is $7^ = (3/2)J76/a, where VL is the Kcplerian orbital frequency 
for a semi-major axis a. If in addition b ^ Rh ~ (m/MQ)^/'^a, the change in eccentricity from 
one encounter is = Ak{m/MQ){b/a)~'^ , where « 6.67 collects the order-unity coefficients 
(Goldreich t Tremaine 1978; Petit & Henon 1986). 

The eccentricity is not the only Kcplerian element that characterizes the non-circular motion 
of a particle; the longitude of periapse specifies the relative orientation of a particle's epicycle. 
The particles may also follow orbits that do not lie in the disk. However, shear-dominated viscous 
stirring excites inclinations at a rate that is always slower than the excitation of eccentricities 
(Wetherill &; Stewart 1993; Goldreich et al. 2004; Rafikov 2003). The perpendicular velocities are, 
in this case, always negligible compared to the epicyclic motion in the disk plane. 

The magnitude of an orbit's eccentricity and the longitude of periapse together specify a 
two-dimensional parameter space. We describe the two-dimensional variable with a vector, e = 
{e cos a;, e sin a;}. The distribution function is a function of this vector and time, f{e,t). That the 
changes in e due to encounters do not depend on the longitude of periapse already shows that the 
distribution function must be axisymmetric, or f{e,t) = f{e,t). Then the number of bodies per 
unit logarithmic interval around e is given by 27re^/(e,t). 

We characterize the eccentricity growth with a differential rate, p{ek)(Pek, that the eccentricity 
vector of a particle will be changed by an amount Cfc. Since the change in eccentricity experienced 

by a pair of bodies, when treated as a vector quantity, is independent of the initial eccentricity 
vector of each body, this function is also axisymmetric and only depends on the magnitude of the 
change of eccentricity, e^. 

The excitation rate depends on the surface mass density of particles in the disk, a, the mass of 
a single body, m, the mass of the central star, AIq, the cross-section at which a particle experiences 
encounters of a strength e^, and the relative speed of those encounters. The impact parameter at 
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which a particle receives an eccentricity scales as 6 ^ . If the eccentricities are small, the 
speed at which one particle encounters the others is set only by the shearing of their two orbits, 
which is proportional to b. Then, as shown in Paper I, 



After simplification, we find 



2'Kp(ek)ekdek = 3^Qb{ek)db{ek). (1) 



Akaa'^ 1 



An integral over every dictates the rate of change of the number of bodies with a given eccen- 
tricity, e: 



df{e,t) 

at 



J Jp{\e-en\)[fien,t)-f{e,t)]d^en. (3) 



Note that this equation implicitly conserves the total particle number, J J f{e,t)d^e = 1. This can 
be shown by integrating both sides with respect to e. 



3. THE SELF- SIMILAR DISTRIBUTION 

Without a specific eccentricity scale to dictate the evolution of f{e,t), we expect a solution of 
the form, 

fie,t) = Fit)g{e/e,it)). (4) 
Replacing f{e,t) in equation 3 with equation 4, we find, 

1 dF{t) 



F{t) at 



^-ec{t)g{x) - ^-^^^c{t) = J J Pi\^ ~ ^n\) [g{xn) - g{x)] c?^x„, (5) 



where x = e/ec{t). The additional constraint that equation 3 conserves particle number implies 
F{t)ec{t)^ is constant. This relationship simplifies the left side of equation 5 such that the only 
possible time-dependence of each term is contained in edt). The right-hand side, however, is 
independent of time. Therefore ec{t) must be constant. Then, 



e^{t) = Cet and F{t) = (Cei)~^ 



(6) 
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The overall normalization of F{t) is arbitrary, as it can be absorbed into g{x). Our choice of 
F{t) requires J J g{x)d?x = 1 to ensure that J J f{e,t)d?e = 1 for all t. Physically, the typical 
eccentricity, edt), is set by the eccentricity change that occurs once per particle per time t, or, 
ecit)'^p{ec{t))t ~ 1. This argument sets edt) only up to a constant coefficient; for simplicity we 
choose the coefficients such that edt) = {Ak/2){aa? /MQ)Q.t. Then, the profile shape, g{x), is 
specified by the integro-differential equation 

Equation 7 is identical to equation 17 of Paper I. A detailed description of the equation and a proof 
of its solution can be found in that paper. For reference, the solution is 

g{x) = ^{l + xY"'\ (8) 



4. NUMERICAL SIMULATIONS 

The non-equilibrium distribution function of eccentricities in the regime discussed above can 
be measured directly from a full numerical simulation of the disk. We use a custom N-body 
integrator that evolves the changes in the two-body constants of motion of each particle around the 
central mass. These constants are chosen to vary slowly with small perturbations. Solving Kepler's 
equation for each body translates each time-step into a change in orbital phase. The constants of 
motion are then integrated by a fourth-order Runge-Kutta routine with adaptive time-steps (Press 
et al. 1992). 

For this study we follow a disk of two-hundred equal mass bodies, with m = 5 x lO^^M©, 
on initially circular orbits with randomly determined phases and semi-major axes within a small 
annulus of width Aa = 0.8a. To avoid possible artifacts from the edge of the simulation, we only 
measure the eccentricities of the bodies in the central third of the disk. A histogram of those 
eccentricities shows the number of bodies with each eccentricity, e dN/de. To increase the signal 
to noise ratio of the histogram at each time, we add the results of one hundred simulations with 
randomly generated initial semi-major axes and orbital phases. 

Figure 1 shows the eccentricity distributions measured after three and ten orbits. The horizon- 
tal error bars indicate the width of each bin, and the vertical error bars are determined assuming 
that each bin is Poisson distributed. The analytic distribution function derived in §3 for each time 
is also plotted, as a solid line. The measured distributions agree remarkably with the analytic 
result. 

To emphasize the self-similarity of the distribution shape, we scale the eccentricities measured 
at each time by the characteristic eccentricity at that time {edt), given by equation 6) and plot the 
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Fig. 1. — Eccentricity distributions of a shear-dominated disk of 200 particles each with a mass 
m = 5 X 1O~^M0 after three (black line) and ten (gray line) orbits. The average surface mass 
density of the simulated annulus is 3 x 10~^g cm~^. The vertical error bars are estimated by 
assuming each bin obeys Poisson statistics. The width of each bin has been chosen such that each 
bin contains a similar number of particles. 
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shapes added together. Figure 2 shows that the resulting distribution shape matches the analytic 
form of g{x) very well. 



We have written a time-dependent Boltzmann equation that describes the eccentricity distri- 
bution function of a population of orbiting particles under the influence of their mutual excitations 
in the shear-dominated regime. Reasoning that the distribution function of eccentricities should 
behave self-similarly, that is, retain a constant profile while its normalization and scaling depend 
on time, we have reduced the Boltzmann equation to a form that has already been shown to have 
an analytic solution. Numerical experiments confirm the self-similarity and the analytic solution. 

Although we have only considered disks of a single particle size, the formalism above applies 
trivially to disks with mass distributions. In fact, the characteristic eccentricity, equation 6, depends 
only on the total surface mass density of the disk. As long as bodies of every part of the mass 
spectrum are in the shear-dominated regime, the eccentricities of all bodies are drawn from the same 
distribution. This is a consequence of the fact that gravitational acceleration is mass- independent. 
In contrast, the dynamical friction of Paper I depends on the size of each particle. The equilibrium 
distributions in that case do differ for each mass group. 

Since ec{t) is an increasing function of time, the condition of shear-dominated dynamics will 
be violated eventually. While the disk is shear-dominated however, most of the disk bodies have 
eccentricities of about edt). The mean eccentricity, / / ef{e,t)(fe, is formally infinite; in reality 
the mean depends logarithmically on the maximum eccentricity achievable from one interaction. 
Higher moments of the distribution, such as (e^), are dominated by the bodies with the maximum 
eccentricity. The random kinetic energy of the disk bodies, for example, is then set by the few 
bodies with the highest eccentricities regardless of the value of edt). 

That the shape of the distribution function is identical to the distribution function of Paper 
I is ultimately not surprising. In both scenarios the bodies in question excite their orbital param- 
eters via the same shear-dominated viscous stirring mechanism. If dynamical friction is acting on 
these bodies, their eccentricities decrease with time proportionally to their magnitude. An equi- 
librium between excitations and this damping produces a characteristic eccentricity around which 
the eccentricities of all bodies are distributed. Without an agent of dynamical friction, the typical 
eccentricity of a body in the disk, etypicai ~ ec(i), grows with time. However, the ratio of the 
eccentricity of a particle that has not interacted recently, e, to that typical eccentricity shrinks 
proportionally to itself: 



This is formally equivalent to the damping provided by the dynamical friction of Paper I. 



5. 



DISCUSSION 




(9) 
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e/e,(t) 

Fig. 2. — The eccentricity distribution of the same numerical simulation of Figure 1. Here, the 
distribution after one, three, and ten orbits, scaled by the characteristic eccentricity at that time 
are added together. The profile shape is very well described by equation 8, plotted as a solid line. 
The error bars are chosen in the same way as Figure 1. 
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